1 Modelo estadístico

El modelo de regresión lineal múltiple generaliza el modelo lineal simple por incluir \(k\) variables regresoras. Se puede formular de dos formas distintas:

\[ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_{p-1} X_{ip-1} + \varepsilon_i,\quad \varepsilon_i \sim N(0, \sigma^2) \label{eq:modelo} \]

donde:

  • \(Y_i\) es la variable respuesta para la observación \(i\), con \(i=1,\ldots,n\).

  • \(\beta_j\) son los coeficientes del modelo, con \(j=0,\ldots,p-1\).

  • \(X_{ij}\) es la \(j\)-ésima covariable para esa observación.

  • \(\varepsilon_i\) es el término de error aleatorio.

  • \(p\) es la cantidad de parámetros a estimar incluyendo el intercept.


1.1 Notación Matricial

En forma matricial, el modelo se escribe como:

\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

donde:

  • \(\mathbf{Y}\) es un vector \(n \times 1\) de respuestas,

  • \(\mathbf{X}\) es una matriz \(n \times p\) de diseño (incluye una columna de 1s),

  • \(\boldsymbol{\beta}\) es un vector \(p \times 1\) de parámetros,

  • \(\boldsymbol{\varepsilon}\) es un vector de errores aleatorios, con \(\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I})\).


Estos vectores y matrices son:

\[\begin{equation} \begin{array}{cc} \underset{n \times 1}{\mathbf{Y}}=\left[\begin{array}{c} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{array}\right] & \underset{n \times p}{\mathbf{X}}=\left[\begin{array}{ccccc} 1 & X_{11} & X_{12} & \cdots & X_{1, p-1} \\ 1 & X_{21} & X_{22} & \cdots & X_{2, p-1} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & X_{n 1} & X_{n 2} & \cdots & X_{n, p-1} \end{array}\right] \\ \underset{p \times 1}{\boldsymbol{\beta}}=\left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{p-1} \end{array}\right] & \underset{n \times 1}{\varepsilon}=\left[\begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{array}\right] \end{array} \end{equation}\]

1.2 Observación

Notemos que, para el caso de dos variables regresoras, el modelo resulta:

\[ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}+\varepsilon_i,\quad con \quad \varepsilon_i \sim N(0, \sigma^2) \]

En este caso, la esperanza condicional de \(Y\) dada las variables \(X_1\) y \(X_2\) es

\[\begin{equation} E(Y|X_{1},X_{2}) = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2}, \label{eq:modelo2} \end{equation}\]

A esta ecuación se la llama superficie de respuesta y su gráfico es un plano. En el caso de regresión lineal simple tenemos que \(E(Y|X) = \beta_0 + \beta_1 X\) y su gráfico es una recta.

2 Supuestos del modelo

  • \(\operatorname{los} \varepsilon_i\) tienen media cero, \(E\left(\varepsilon_i\right)=0\).

  • \(\operatorname{los} \varepsilon_i\) tienen todos la misma varianza desconocida que llamaremos \(\sigma^2\) y que es el otro parámetro del modelo, \(\operatorname{Var}\left(\varepsilon_i\right)=\sigma^2\).

  • \(\operatorname{los} \varepsilon_i\) tienen distribución normal.

  • \(\operatorname{los} \varepsilon_i\) son independientes entre sí, e independientes de las covariables \(X_{i 1}, X_{i 2}\), \(\ldots, X_{i p-1}\).

3 Estimación de los parámetros del modelo

3.1 Estimador de mínimos cuadrados

\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \]

  • Este estimador minimiza la suma de cuadrados de los residuos dada por

\[ \min_{\boldsymbol{\beta}} \| \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \|^2 \]

3.2 Estimador por máxima verosimilitud

De la misma manera que sucede en regresión lineal simple, se puede probar que el estimador de máxima verosimilitud para los parámetros del modelo coinciden con los de mínimos cuadrados, siempre que los errores se distribuyan en forma normal y sean independientes.

  • El estimador de \(\sigma^2\) es

\[ \hat{\sigma}^2=\dfrac{\sum_{i=1}^n (Y_{i}-\widehat{Y}_{i})^2}{n-p}\] donde \(p\) es la cantidad de parámetros a estimar.

4 Ejemplo

Un embotellador de bebidas gaseosas quiere mejorar la eficiencia de sus rutas de distribución. Para eso, está estudiando cuánto tiempo le toma a un representante de ruta atender las máquinas expendedoras en una tienda. Este proceso de atención incluye reponer productos embotellados en la máquina y realizar tareas básicas de limpieza y mantenimiento. Un ingeniero industrial que colabora en el estudio propone que el tiempo que se tarda en atender una máquina depende principalmente de dos factores: la cantidad de cajas de producto que se reponen (\(X_1\)) y la distancia que el representante tiene que caminar dentro de la tienda (\(X_2\) medida en metros). Se han recogido 25 observaciones sobre el tiempo de servicio.

Este ejemplo se obtuvo del libro Introducción al análisis de regresión lineal de Montgmoery et al. que se encuentra en las referencias.

4.1 Lectura del conjunto de datos

Tiempos <- read.csv("./Datos/Tiempos1.csv")

4.2 Exploración de las relaciones entre las variables

library(ggplot2)
library(GGally)

ggpairs(Tiempos,
        columns = 1:3,
        mapping = aes(alpha = 0.9),
        upper = list(continuous = wrap("points", color = "darkred")),
        lower = list(continuous = wrap("points", color = "blue")),
        diag = list(continuous = wrap("densityDiag", fill = "lightblue")))

4.3 Scatter plot en 3D

library(scatterplot3d)
scatterplot3d(x=Tiempos$cajas, y=Tiempos$distancia, z=Tiempos$tiempo, pch=16, cex.lab=1,
              highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",
              ylab="Distancia", zlab="Tiempo")

Se puede observar que, a medida que la cantidad de cajas y la distancia aumentan, los tiempos también aumentan.

Se puede realizar un diagrama de dispersión con movimiento.

library(plotly)
plot_ly(x=Tiempos$cajas, y=Tiempos$distancia, z=Tiempos$tiempo, type="scatter3d", color=Tiempos$tiempo) %>% 
  layout(scene = list(xaxis = list(title = "Cantidad"),
                      yaxis = list(title = "Distancia"),
                      zaxis = list(title = "Tiempo")))

4.4 El modelo

\[ Y_i=\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\varepsilon_i,\quad con \quad \varepsilon_i \sim N(0, \sigma^2) \] donde

  • \(Y\): tiempo de atención a las máquinas expendedoras.

  • \(X_1\): cantidad de cajas entregadas.

  • \(X_2\): distancia recorrida por el representante.

4.4.1 El modelo en R

modelo <- lm(tiempo ~ cajas + distancia, data=Tiempos)

4.4.2 La salida del modelo

summary(modelo)
## 
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9508 -1.7810  0.0006  1.8293  4.3545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.137487   0.986207   5.209 3.66e-05 ***
## cajas       1.526082   0.134819  11.320 2.12e-10 ***
## distancia   0.008508   0.002955   2.879  0.00897 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.516 on 21 degrees of freedom
## Multiple R-squared:  0.943,  Adjusted R-squared:  0.9376 
## F-statistic: 173.7 on 2 and 21 DF,  p-value: 8.666e-14

4.4.3 Modelo propuesto

\[ E(Y|X_{1},X_{2}) =2.341231+1.615907X_1+0.014385 X_2 \]

4.4.4 El scatter plot junto con la superficie de respuesta

# Se crea el grafico 3d y se guarda en un objeto, por ejemplo mi_3d
mi_3d <- scatterplot3d(x=Tiempos$cajas, y=Tiempos$distancia, z=Tiempos$tiempo, pch=16,            cex.lab=1,highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",ylab="Distancia ", zlab="Tiempo")

# Para agregar el plano usamos $plane3d( ) con argumento modelo ajustado
mi_3d$plane3d(modelo, lty.box="solid", col="mediumblue")

4.5 Interpretación de los coeficientes del modelo

  • \(\beta_0=2.341231\) indica el tiempo promedio estimado de servicio cuando:

    • no se entregan cajas \(X_1=0\),

    • no se camina distancia alguna \(X_2=0\).

  • \(\beta_1 = 1.615907\) indica que el tiempo de servicio aumenta, en promedio, aproximadamente \(2.34\) minutos cuando se agrega una caja adicional y se mantiene la distancia constante.

  • \(\beta_2 = 0.014385\) indica que el tiempo de servicio aumenta, en promedio, aproximadamente \(0.01\) minutos cuando la distancia recorrida aumenta un metro y se mantiene la cantidad de cajas constante.

5 Validación de los supuestos

5.1 Residuos vs. valores ajustados

ggplot(data = Tiempos, aes(x = fitted(modelo), y = rstandard(modelo))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "valores ajustados", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Valores predichos")

5.2 Residuos vs. cada variable regresora

Estos gráficos sirven para detectar la existencia de una relación curvilínea entre los residuos y cada variable regresora. Esto indicaría que debiera incluirse un término \(X_i^2\) o \(ln(X_i)\) o alguna otra transformación.

ggplot(data = Tiempos, aes(x = cajas, y = rstandard(modelo))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "cantidad de cajas", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Cantidad de cajas")

ggplot(data = Tiempos, aes(x = distancia, y = rstandard(modelo))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "distancia", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Distancia")

5.3 Normalidad de los residuos estandarizados

  • Qqnorm
residuos.df=data.frame(resid=rstandard(modelo))
ggplot(data = residuos.df, aes(sample = rstandard(modelo))) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "QQ Plot de los residuos",
       x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
  theme_minimal()

  • Boxplot
ggplot(residuos.df, aes(y = rstandard(modelo))) +
  geom_boxplot(fill = "skyblue") +
  labs(title = "Boxplot de residuos estandarizados",
       y = "Residuos estandarizados") +
  theme_minimal()

6 Medidas de bondad del ajuste del modelo: \(R^2\) y \(R^2_{a}\)

6.1 \(R^2\)

  • Recordemos la definición de \(R^2\), la medida que indica la proporción de variablidad explicada por el modelo.

\[ R^2=1-\frac{\text { Suma de cuadrados del residuo ($SSRes$) }}{\text { Suma total de cuadrados (SST) }} \] - Este coeficiente tiene la misma interpretación que en regresión lineal simple.

  • Pero…, hay que notar que si agregamos variables al modelo, el valor de \(SSRes\) disminuye.

  • Es decir, si agregamos variables al modelo, el valor de \(R^2\) aumenta en forma espúrea.

  • Por este motivo,cuando queremos comparar modelos de regresión que tienen distinta cantidad de variables, no conviene usar el \(R^2\),

  • Conviene usar otra medida que corrige ese valor teniendo en cuenta la cantidad de variables que se usaron en el modelo.

6.2 \(R^2\) ajustado (\(R^2_a\))

El \(R^2_a\) penaliza por el número de variables en el modelo y solo aumenta si la nueva variable mejora realmente el modelo. Se define:

\[ R_{a}^2=1-\left(\frac{\left(1-R^2\right)(n-1)}{n-k-1}\right) \]

donde:

  • \(n\) es el número de observaciones,

  • \(k\) es el número de predictores.

7 Test de hipótesis para los parámetros del modelo

Se puede analizar si cada parámetro es indivdualmente significativamente distinto de cero, de la misma forma que hicimos en regresión lineal simple. Los test \(\textit{t}\) sirven para contrastar si cada parámetro individual \(\textit{\beta_j}\) es significativamente distinto de 0. Es decir, testear las siguientes hipótesis:

\[ H_0: \ \beta_j=0 \text{ vs. } \ H_1: \ \beta_j \neq 0 \ \text{ con } \ j=1,2\] Donde:

  • \(\beta_1\) es el efecto de la variable \(\textit{cajas}\) (que mide el número de productos entregados).

  • \(\beta_2\) es el efecto de la variable \(\textit{distancia}\) (que mide la distancia caminada).

7.1 Estadísticos de contraste

\[T=\dfrac{\hat{\beta}_j}{se(\hat{\beta}_j)} \sim t_{n-p-1} \ \text{ bajo } \ H_0\]

  • Se rechaza \(H_0\) si \(|t_c|>t_{\alpha/2,n-p-1}.\)

  • Es importante notar que esta es una prueba parcial o condicional, ya que el coeficiente de regresión \(\beta_j\) depende de la presencia de todas las demás variables explicativas \(x_i\) (con \(i \neq j\)) en el modelo. Por lo tanto, el test \(t\) para \(\beta_j\) evalúa el aporte de la variable \(x_j\) al modelo, condicional a que las demás variables ya están incluidas. En otras palabras, nos dice si \(x_j\) mejora significativamente la explicación de la variable respuesta, dado que las demás regresoras ya están en el modelo.

7.2 Interpretación de la salida summary(modelo):

summary(modelo)
## 
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9508 -1.7810  0.0006  1.8293  4.3545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.137487   0.986207   5.209 3.66e-05 ***
## cajas       1.526082   0.134819  11.320 2.12e-10 ***
## distancia   0.008508   0.002955   2.879  0.00897 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.516 on 21 degrees of freedom
## Multiple R-squared:  0.943,  Adjusted R-squared:  0.9376 
## F-statistic: 173.7 on 2 and 21 DF,  p-value: 8.666e-14
  • Todos los p-valores son chicos, quiere decir que se rechaza la hipótesis nula a cualquier nivel indicando que cada variable aporta información explicativa más allá del azar.

  • Por lo tanto, tanto la variable \(\textit{cajas}\) como la variable \(\textit{distancia}\) influyen significativamente en el tiempo de entrega.

  • Si algún p-valor es grande, no se rechaza la hipótesis nula indicando que la/s variable/s no influyen significativamente en el modelo, pudiendo también indicar presencia de colinealidad (ver VIF más abajo).

8 Prueba de significancia de la regresión - Test F

Este test permite evaluar si el modelo es globalmente significativo, es decir, si al menos una de las variables regresoras explica la variabilidad de la respuesta.

  • Este test se basa en el cociente de la mejora debida al modelo (SSReg) y la diferencia entre el modelo y los datos observados (SSRes).

  • El objetivo es ver si al menos uno de los coeficientes \(\beta_j\) de los predictores es significativamente distinto de cero.

8.1 Hipótesis del test F

\[H_0: \beta_1=\ldots=\beta_{p}=0 \text{ vs. } H_1: \beta_j=0 \text{ para algún } j=1,\ldots,p\]

  • Si se rechaza \(H_0\) implica que al menos una de las variables regresoras contribuye al modelo en forma significativa.

8.2 Expresión del estadístico F

  • Recordemos

\[SS_T=SSReg+SSRes \]

  • Entonces

\[\text{F}=\dfrac{\dfrac{\text{Variación explicada por el modelo}}{p-1}}{\dfrac{\text{Variación residual}}{n−p}}=\dfrac{{\dfrac{SSReg}{p-1}}}{{\dfrac{SSRes}{n-p}}}=\dfrac{MSReg}{MSRes}.\]

Donde:

  • \(n\): número de observaciones
  • \(p\): número de variables predictoras
  • \(SSReg\): suma de cuadrados de la regresión
  • \(SSRes\): suma de cuadrados del error

8.3 Distribución del estadístico F.

\[\text{F}\sim F_{p-1,n-p}\]

  • Regla de decisión:

    • Se rechaza \(H_0\) a nivel \(\alpha\) cuando \(\text{F} >F_{p-1,n-p,\alpha}\), donde \(F_{p-1,n-p,\alpha}\) es el valor del estadístico que deja un área de \(\alpha\) a la derecha de él, cuando \(H_0\) es verdadera.
  • Valores grandes de F corresponde a p-valores pequeños.

  • Si se rechaza \(H_0\) indica que no todos los coeficientes del modelo son nulos.

  • Un valor alto de \(F\) indica que el modelo es significativo.

8.3.1 Ejemplo

summary(modelo)
## 
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9508 -1.7810  0.0006  1.8293  4.3545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.137487   0.986207   5.209 3.66e-05 ***
## cajas       1.526082   0.134819  11.320 2.12e-10 ***
## distancia   0.008508   0.002955   2.879  0.00897 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.516 on 21 degrees of freedom
## Multiple R-squared:  0.943,  Adjusted R-squared:  0.9376 
## F-statistic: 173.7 on 2 and 21 DF,  p-value: 8.666e-14
  • En el último renglón del summary vemos que

    • F-statistic: 173.7 on 2 and 21 DF, p-value: 8.666e-14.
  • Por lo tanto el pvalor es muy chico, se rechaza \(H_0\) a nivel \(5\%\) indicando que, al menos algún \(\beta_j\) es no nulo.

  • Recordemos

\[ p-\text { valor }=P\left(F_{p-1, n-p}>F_{obs}\right) \]

p=3
n=nrow(Tiempos)
n
## [1] 24
pf(173.7,p-1, n-p,lower.tail=FALSE)
## [1] 8.648615e-14

8.3.2 ¿Qué son los grados de libretad?

  • al ajustar el modelo estamos usando \(p\) valores (los coeficientes beta) que se determinan a partir de los datos.

  • esos valores se obtienen resolviendo \(p\) ecuaciones, que se conocen como las ecuaciones normales.

  • estas ecuaciones surgen al minimizar el error del modelo: se derivan \(p\) veces (una por cada parámetro) y se igualan a cero.

  • entonces, los residuos no son completamente libres: están “condicionados” por esas \(p\) ecuaciones.

  • luego, si conocemos \(n-p\) residuos, los otros \(p\) se pueden calcular a partir de ellos.

  • por eso, decimos que los residuos tienen \(n-p\) grados de libertad.

En un modelo de regresión lineal múltiple, los coeficientes estimados representan el cambio esperado en la variable respuesta asociado con una unidad de cambio en cada predictor, manteniendo los demás constantes. Sin embargo, estas estimaciones son sujetas a variabilidad muestral.

Para cuantificar la precisión de estos estimadores, se calculan intervalos de confianza (IC), que representan rangos plausibles donde se espera que se encuentre el valor real del parámetro con un cierto nivel de confianza (usualmente 95%).

Matemáticamente, un intervalo de confianza al 95% para un coeficiente \(\beta_j\) se calcula como:

\[ IC_{95\%}(\beta_j) = \hat{\beta}_j \pm t^* \cdot SE(\hat{\beta}_j) \]

donde: - \(\hat{\beta}_j\) es la estimación puntual del coeficiente, - \(SE(\hat{\beta}_j)\) es el error estándar de dicha estimación, - \(t^*\) es el valor crítico de la distribución t de Student con \(n - p - 1\) grados de libertad (número de observaciones menos número de parámetros).

Un intervalo que incluye el cero sugiere que no hay evidencia suficiente para afirmar que el coeficiente sea distinto de cero al nivel de confianza elegido.

confint(modelo)
##                   2.5 %     97.5 %
## (Intercept) 3.086557749 7.18841554
## cajas       1.245711004 1.80645228
## distancia   0.002363293 0.01465241

9 Análisis de puntos influyentes

9.1 Outliers

Es importante destacar que, aunque las observaciones \((X_{i1}, X_{i2}, \dots, X_{i(p-1)}, Y_i)\) no pueden representarse gráficamente en más de dos dimensiones, siempre es posible analizar gráficamente:

  • Los valores ajustados \(\hat{Y}_i\)
  • Los residuos \(r_i\)

A continuación, se listan gráficos útiles para evaluar la idoneidad del modelo de regresión múltiple:

9.1.1 Gráfico de residuos vs. valores ajustados

  • Permite evaluar:
    • La adecuación del modelo a los datos observados
    • La homoscedasticidad (detectar errores con varianza no constante)
    • La presencia de valores extremos (outliers)

9.1.2 Gráfico de residuos vs. tiempo (u orden de recolección de datos)

  • Útil cuando los datos tienen un orden temporal o secuencial
  • Permite detectar correlación entre los errores (autocorrelación)

9.1.3 Boxplots y gráfico de probabilidad normal de los residuos

  • Sirven para evaluar si los errores:
    • Se distribuyen normalmente
    • Tienen simetría o presentan colas pesadas

9.1.4 Gráficos de residuos vs. cada variable predictora

  • Permiten evaluar:
    • Si la relación entre la variable predictora y la respuesta es lineal
    • Si se podría requer un término cuadrático o alguna transformación
    • Si la varianza del error depende de dicha variable

9.1.5 Gráficos de residuos vs. variables omitidas importantes

  • Ayudan a detectar efectos importantes no incluidos en el modelo
  • Pueden revelar relaciones adicionales con la variable de respuesta

9.1.6 Gráficos de residuos vs. términos de interacción

  • Evaluación visual para decidir si se necesitan términos de interacción
  • Estos términos pueden capturar efectos combinados entre predictores

9.2 Leverage

Recordemos que el leverage mide qué tan lejos están las observaciones respecto del centro de los datos.

  • ¿Cuándo una observación \(i\) se considerará que tiene alto leverage? Cuando se cumple que

\[h_{ii}>2\bar{h}\]

  • En regresión lineal simple, \(\sum_{i=1}^n h_{ii}=2\). Por lo tanto, \(\bar{h}=\dfrac{2}{n}.\) Entonces, una observación se considerará que tiene alto leverage si \[h_{ii}>\dfrac{4}{n}.\]

  • En regresión lineal múltiple, \(\sum_{i=1}^n h_{ii}=p\). Por lo tanto, \(\bar{h}=\dfrac{p}{n}.\) Entonces, una observación se considerará que tiene alto leverage si \[h_{ii}>\dfrac{2 p}{n}\]

  • Cabe señalar, que este criterio está ajustado por la cantidad de covariables.

  • Algunos autores utilizan otro criterio sin considerar este ajuste.

    • Se considera que una observación tiene alto leverage si \[h_{ii}>0.5.\]

    • Se considera que una observación tiene un moderado leverage si \[0.2<h_{ii}<0.5.\]

  • También se puede graficar un histograma de los valores \(h_{ii}\) y observar si hay algún valor/es que se comporta/n encuentran lejanos al resto.

9.2.1 ¿Cómo detectarlos en R?

hatvalues(modelo)
##          1          2          3          4          5          6          7 
## 0.11992977 0.10097319 0.10024969 0.10128434 0.16642759 0.08066357 0.07620000 
##          8          9         10         11         12         13         14 
## 0.06254612 0.06537738 0.11367920 0.07741039 0.06437975 0.04290146 0.09998709 
##         15         16         17         18         19         20         21 
## 0.04661503 0.04687996 0.08033747 0.11083391 0.20438000 0.13894064 0.14675966 
##         22         23         24 
## 0.18537865 0.21115081 0.55671434
p=3
which(hatvalues(modelo) > 2*p/dim(Tiempos)[1])
## 24 
## 24

9.2.2 Histograma de \(h_{ii}\)

leverage.df<-data.frame(x=hatvalues(modelo))

ggplot(leverage.df, aes(x=x))+
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

9.3 Cook’s distance

  • Mide la influencia que tiene una observación determinada en la estimación de los parámetros de la regresión.

  • Una observación con distancia de Cook grande indica que la ausencia de esa observación podría cambiar los valores estimados de los parámetros.

  • Es una medida que toma en cuenta tanto los residuos como los puntos de alto leverage.

\[ D_i = \frac{\displaystyle{\sum_{j=1}^n}\left(\widehat{Y}_j - \widehat{Y}_{(i)j}\right)^2}{p \widehat{\sigma}^2} \] donde: - \(p\) es la cantidad de parámetros a estimar en el modelo.

  • \(\widehat{Y}_{(i)j}\) corresponde al valor predicho para la j-ésima observación si se usaron las \(n-1\) restantes observaciones para hacer el ajuste.

  • \(\widehat{Y}_{j}\) es el valor predicho para la j-ésima observación en el modelo ajustado con las \(n\) observaciones.

9.4 Otra forma de calcular los \(D_i\).

\[ D_i = \dfrac{r_{si} ^2}{p} \frac{h_{ii}}{1 - h_{ii}} \] donde \(r_{si}\) es el i-ésimo residuo standarizado y \(h_{ii}\) es el i-ésimo valor del leverage, ambos definidos anteriormente.

El \(r_{si}\) mide cuán lejos está el valor predicho de la observación, mientras que el \(h_{ii}\) mide el grado de apalancamiento de la observación \(i\). Así, un valor elevado de \(D_{i}\) puede deberse a un gran valor de \(r_{i}\) , a un gran valor de \(h_{ii}\) o a ambos.

Entonces:

  • \(D_i\) depende de dos factores:

    • los \(r_i\) y los \(h_{ii}\).

    • Cuánto más grande son los \(r_i\) y los \(h_{ii}\), mayor será el valor de \(D_i\).

  • Por lo tanto, la \(j-ésima\) observación puede ser influyente si sucede alguna de las siguientes situaciones:

    • tener un \(r_i\) alto y sólo un moderado valor de \(h_{ii}\).

    • tener un valor alto valor de \(h_{ii}\) con sólo un moderado valor de \(r_i\).

    • tener tanto un alto valor de \(r_i\) como un alto valor de \(r_i\).

9.5 Criterio para determinar si una observación tiene un \(D_i\) alto.

Los valores de \(D_i\) se comparan con los percentiles de la distribución \(F\) de Fisher con \(p\) y \(n-p\) grados de libertad en el numerador y denominador, respectivamente.

El criterio para decidir si una observación es influyente es el siguiente:

  • Si \(D_i<\) percentil 0,20 de la distribución \(F_{p, n-p}\) entonces la observación no es influyente.

  • Si \(D_i>\) percentil 0,50 de la distribución \(F_{p, n-p}\) entonces la observación es muy influyente y requerirá tomar alguna medida.

  • Si \(D_i\) se encuentra entre el percentil 0,20 y el percentil 0,50 de la distribución \(F_{p, n-p}\) se sugiere mirar además otros estadisticos.

9.5.1 ¿Cómo detectarlos en R?

En el ejemplo:

n=dim(Tiempos)[1]
p=3
limite05.Cook=qf(0.5,p, n-p)
limite05.Cook
## [1] 0.8148716
limite02.Cook=qf(0.2,p, n-p)
limite02.Cook
## [1] 0.3352225
cooks.distance(modelo)
##            1            2            3            4            5            6 
## 2.221728e-03 5.753886e-02 8.259894e-03 2.793573e-02 1.716661e-01 1.664864e-02 
##            7            8            9           10           11           12 
## 1.469579e-02 7.105612e-02 1.964887e-04 2.393920e-02 3.662926e-05 3.045623e-05 
##           13           14           15           16           17           18 
## 6.607918e-04 2.174400e-03 4.678501e-03 4.241174e-02 1.163750e-02 1.125283e-01 
##           19           20           21           22           23           24 
## 2.184760e-01 8.405315e-02 9.061495e-02 9.429366e-02 7.143746e-02 5.606107e-02

Hay varias formas.

# Cook's Distance
plot(modelo, which = 4)

# Distancia de Cook
cooksd <- cooks.distance(modelo)

# Ver los índices de los puntos más influyentes
which(cooksd >limite05.Cook)
## named integer(0)
which(limite02.Cook<cooksd)
## named integer(0)

10 Análisis de las observación influyentes

  • Defino el conjunto de datos sin la observación \(24\).
Tiempos.sin<-Tiempos[-24,]
  • Defino el modelo
modelo.sin<-lm(tiempo ~ cajas + distancia, data=Tiempos.sin)
  • Salida del modelo
summary(modelo.sin)
## 
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos.sin)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9335 -1.8143 -0.2151  2.0083  4.3673 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.347585   1.165595   4.588 0.000178 ***
## cajas       1.480844   0.186832   7.926 1.35e-07 ***
## distancia   0.008750   0.003093   2.829 0.010364 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.57 on 20 degrees of freedom
## Multiple R-squared:  0.8951, Adjusted R-squared:  0.8846 
## F-statistic: 85.35 on 2 and 20 DF,  p-value: 1.61e-10
  • Gráfico de Residuos vs. valores ajustados
ggplot(data = Tiempos.sin, aes(x = fitted(modelo.sin), y = rstandard(modelo.sin))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "valores ajustados", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Valores predichos")

  • Gráfico de Residuos vs. cada variable regresora
ggplot(data = Tiempos.sin, aes(x = cajas, y = rstandard(modelo.sin))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "cantidad de cajas", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Cantidad de cajas")

ggplot(data = Tiempos.sin, aes(x = distancia, y = rstandard(modelo.sin))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "distancia", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Distancia")

  • Normalidad de los residuos

    • Qqnorm
residuos.sin.df=data.frame(resid=rstandard(modelo.sin))

ggplot(data = residuos.sin.df, aes(sample = rstandard(modelo.sin))) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "QQ Plot de los residuos",
       x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
  theme_minimal()

  • Boxplot
ggplot(residuos.sin.df, aes(y = rstandard(modelo.sin))) +
  geom_boxplot(fill = "skyblue") +
  labs(title = "Boxplot de residuos estandarizados",
       y = "Residuos estandarizados") +
  theme_minimal()

10.0.1 Conclusión

No se detectan cambios grande en la estimación de los parámetros, la estimación de la varianza de los residuos, en las decisiones de los test individuales ni en el test global. Por lo tanto se elige el modelo con la observación \(24\).

11 Otro ejemplo

11.1 Lectura del conjunto de datos

Tiempos2 <- read.csv("./Datos/Tiempos.csv")

11.2 Exploración de las relaciones entre las variables

library(ggplot2)
library(GGally)

ggpairs(Tiempos2,
        columns = 1:3,
        mapping = aes(alpha = 0.9),
        upper = list(continuous = wrap("points", color = "darkred")),
        lower = list(continuous = wrap("points", color = "blue")),
        diag = list(continuous = wrap("densityDiag", fill = "lightblue")))

11.3 Scatter plot en 3D

library(scatterplot3d)
scatterplot3d(x=Tiempos2$cajas, y=Tiempos2$distancia, z=Tiempos2$tiempo, pch=16, cex.lab=1,
              highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",
              ylab="Distancia", zlab="Tiempo")

12 El modelo en R

modelo2 <- lm(tiempo ~ cajas + distancia, data=Tiempos2)

12.1 La salida del modelo

summary(modelo2)
## 
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7880 -0.6629  0.4364  1.1566  7.4197 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.341231   1.096730   2.135 0.044170 *  
## cajas       1.615907   0.170735   9.464 3.25e-09 ***
## distancia   0.014385   0.003613   3.981 0.000631 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.259 on 22 degrees of freedom
## Multiple R-squared:  0.9596, Adjusted R-squared:  0.9559 
## F-statistic: 261.2 on 2 and 22 DF,  p-value: 4.687e-16

13 El scatter plot junto con la superficie de respuesta

# Se crea el grafico 3d y se guarda en un objeto, por ejemplo mi_3d
mi_3d <- scatterplot3d(x=Tiempos2$cajas, y=Tiempos2$distancia, z=Tiempos2$tiempo, pch=16,            cex.lab=1,highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",ylab="Distancia ", zlab="Tiempo")

# Para agregar el plano usamos $plane3d( ) con argumento modelo ajustado
mi_3d$plane3d(modelo2, lty.box="solid", col="mediumblue")

  • Se puede ver la existencia de un punto influyente.

14 Validación de los supuestos

14.1 Residuos vs. valores ajustados

ggplot(data = Tiempos2, aes(x = fitted(modelo2), y = rstandard(modelo2))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "valores ajustados", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Valores predichos")

14.2 Residuos vs. cada variable regresora

Estos gráficos sirven para detectar la existencia de una relación curvilínea entre los residuos y cada variable regresora. Esto indicaría que debiera incluirse un término \(X_i^2\) o \(ln(X_i)\) o alguna otra transformación.

ggplot(data = Tiempos2, aes(x = cajas, y = rstandard(modelo2))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "cantidad de cajas", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Cantidad de cajas")

ggplot(data = Tiempos2, aes(x = distancia, y = rstandard(modelo2))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "distancia", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Distancia")

14.3 Normalidad de los residuos

residuos.df=data.frame(resid=rstandard(modelo2))
ggplot(data = residuos.df, aes(sample = rstandard(modelo2))) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "QQ Plot de los residuos",
       x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
  theme_minimal()

14.4 Análisis de puntos influyentes

14.5 Outliers

sort(rstandard(modelo2))
##          20           1          24          22          23          21 
## -1.87354643 -1.62767993 -1.49605875 -1.44999541 -1.44368977 -0.87784258 
##          16          12           5           6          25           3 
## -0.22270023 -0.19325733 -0.14176094 -0.09080847 -0.06750861 -0.01609165 
##          17          15           7          13          14           2 
##  0.13803929  0.21029137  0.27042496  0.32517935  0.34113547  0.36484267 
##           8          19          11          10          18           4 
##  0.36672118  0.57876634  0.71807970  0.81325432  1.11295196  1.57972040 
##           9 
##  3.21376278

14.5.1 Leverage

  • Encuentro el leverage de cada observación
leverage <- hatvalues(modelo)
leverage
##          1          2          3          4          5          6          7 
## 0.11992977 0.10097319 0.10024969 0.10128434 0.16642759 0.08066357 0.07620000 
##          8          9         10         11         12         13         14 
## 0.06254612 0.06537738 0.11367920 0.07741039 0.06437975 0.04290146 0.09998709 
##         15         16         17         18         19         20         21 
## 0.04661503 0.04687996 0.08033747 0.11083391 0.20438000 0.13894064 0.14675966 
##         22         23         24 
## 0.18537865 0.21115081 0.55671434
  • Defino el umbral
p<-3
n=nrow(Tiempos2)
leverage_umbral<- 2 * p / n
  • Encuentro las observaciones que superan el umbral
leverage[leverage>leverage_umbral]
##        24 
## 0.5567143

14.6 Distancia de Cook

  • Encuentro las distancia de Cook de cada observación
cooks <- cooks.distance(modelo)
  • Defino el umbral
p <- length(coef(modelo)) 
n=nrow(Tiempos2)
cooks_umbral <- qf(0.5,p,n-p)
  • Encuentro las observaciones que superan el umbral
cooks[cooks>cooks_umbral]
## named numeric(0)
  • \(\textit{Observación:}\) Las observaciones \(9\) y \(24\) son las que presentan residuos altos y alto leverage, respectivamente. Sin embargo, no presentan una distancia de Cook detectable como punto influeyente. Es de buena práctica estudiar el modelo sin estas observaciones para analizar los cambios que se pudieran producir.

15 Multicolinealidad

15.1 ¿Qué es la multicolinealidad?

La multicolinealidad ocurre cuando existe correlación entre dos o más variables regresoras.

15.2 Causas de la multicolinealidad

  • Relaciones lógicas o matemáticas entre variables
    Variables que están construidas a partir de otras.
    Ejemplo: población total y población urbana.

  • Variables derivadas unas de otras
    Incluir tanto una variable como una versión transformada o combinada.
    Ejemplo: incluir tanto peso como el índice de masa corporal (IMC).

  • Variables naturalmente correlacionadas por contexto
    Algunas variables tienden a moverse juntas por razones sociales o económicas.
    Ejemplo: nivel educativo y salario.

  • Redundancia por variables similares
    Incluir indicadores que miden aspectos muy parecidos.
    Ejemplo: gasto mensual y gasto anual.

  • Tamaño de muestra pequeño
    Con pocos datos, las correlaciones entre variables explicativas se amplifican, generando efectos de colinealidad incluso con relaciones moderadas.

15.3 Problemas que causa la multicolinealidad

  • Los estimadores de los coeficientes pueden cambiar drásticamente con pequeños cambios en los datos.

  • Cuando se incluyen covariables altamente correlacionadas en un modelo, los errores estándar de los coeficientes estimados tienden a aumentar de manera artificial. A este fenómeno se lo conoce como inflación de la varianza de los estimadores.

  • Los coeficientes pueden no resultar estadísticamente significativos, incluso cuando existe una relación real entre la variable dependiente y las variables explicativas.

15.4 Medida de diagnóstico de multicolinealidad

  • Un método para detectar colinealidad entre las variables regresoras es el Factor de Infalción de la Varianza (VIF por sus siglas en inglés).

  • Este factor se calcula para cada variable regresora.

  • El VIF para la k-ésima variable regresora se define:

\[\mathrm{VIF}_k=\dfrac{1}{1-R_k^2} \] con \(k=1,\ldots,p-1\) y \(R_k^2\) es el coeficiente de determinación múltiple cuando \(X_k\) es regresado en las \(p- 2\) restantes covariables \(X\) en el modelo.

  • \(\mathrm{VIF}_k=1\) cuando \(R_k^2=0\) : la covariable \(X_k\) no está correlacionada con las demás.

  • Si \(R_k^2>0\), entonces \(\mathrm{VIF}_k>1\) : hay algún grado de colinealidad.

  • Si \(R_k^2 \approx 1\), el \(\mathrm{VIF}_k\) puede ser muy grande, indicando alta colinealidad.

  • \(\mathrm{VIF} > \mathbf{10} \rightarrow\) indicio fuerte de multicolinealidad.

  • \(\mathbf{5}<\mathrm{VIF} < \mathbf{10}\) indicio de moderada multicolinealidad.

  • Otro criterio: si el promedio de los VIF es significativamente mayor a 1, también sugiere problemas de colinealidad.

15.5 Ejemplo

A continuación se muestra un ejemplo comparando dos modelos de regresión:

  • Uno sin multicolinealidad, con tres variables independientes no correlacionadas.

  • Otro con multicolinealidad, donde una variable es casi una copia de otra.

set.seed(123)

# N° de observaciones
n <- 100

# Modelo sin multicolinealidad
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y1 <- 3 + 2 * x1 + 1.5 * x2 + 1 * x3 + rnorm(n)

datos1 <- data.frame(y = y1, x1 = x1, x2 = x2, x3 = x3)


# Modelo con multicolinealidad (x3 ≈ x1)
x3c <- x1 + rnorm(n, 0, 0.05)  # Alta correlación con x1
y2 <- 3 + 2 * x1 + 1.5 * x2 + 1 * x3c + rnorm(n)

y2 <- 3 + 2 * x1 + 1.5 * x2 + 1 * x3c + rnorm(n)

datos2 <- data.frame(y = y2, x1 = x1, x2 = x2, x3 = x3c)

15.6 Modelo sin multicolinealidad

  • Salida del summary
modelo_sin_col <- lm(y ~ x1 + x2 + x3, data = datos1)
summary(modelo_sin_col)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.49138 -0.65392  0.05664  0.67033  2.53210 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.9807     0.1073  27.768  < 2e-16 ***
## x1            1.9445     0.1169  16.638  < 2e-16 ***
## x2            1.5462     0.1095  14.126  < 2e-16 ***
## x3            0.9426     0.1122   8.399 4.04e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.052 on 96 degrees of freedom
## Multiple R-squared:  0.8392, Adjusted R-squared:  0.8342 
## F-statistic:   167 on 3 and 96 DF,  p-value: < 2.2e-16
  • Residuos vs. valores ajustados
ggplot(data = datos1 , aes(x = fitted(modelo_sin_col), y = rstandard(modelo_sin_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "valores ajustados", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Valores predichos")

  • Residuos vs. cada variable regresora
ggplot(data = datos1 , aes(x = x1, y = rstandard(modelo_sin_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = expression(x[1]), y = "residuos estandarizados",
       title = expression(paste("Residuos estandarizados vs. ", x[1])))

ggplot(data = datos1 , aes(x = x2, y = rstandard(modelo_sin_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = expression(x[2]), y = "residuos estandarizados",
       title = expression(paste("Residuos estandarizados vs. ", x[2])))

ggplot(data = datos1 , aes(x = x3, y = rstandard(modelo_sin_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = expression(x[3]), y = "residuos estandarizados",
       title = expression(paste("Residuos estandarizados vs. ", x[3])))

15.7 Normalidad de los residuos

residuos.df=data.frame(resid=rstandard(modelo_sin_col))
ggplot(data = residuos.df, aes(sample = rstandard(modelo_sin_col))) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "QQ Plot de los residuos",
       x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
  theme_minimal()

15.8 Modelo con multicolinealidad

  • Summary
modelo_con_col <- lm(y ~ x1 + x2 + x3, data = datos2)
summary(modelo_con_col)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = datos2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.43303 -0.74973  0.04546  0.70439  2.40904 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.8625     0.1058  27.047   <2e-16 ***
## x1            2.0176     2.1563   0.936    0.352    
## x2            1.5811     0.1094  14.452   <2e-16 ***
## x3            0.9474     2.1768   0.435    0.664    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.04 on 96 degrees of freedom
## Multiple R-squared:  0.8974, Adjusted R-squared:  0.8942 
## F-statistic:   280 on 3 and 96 DF,  p-value: < 2.2e-16
  • Residuos vs. valores ajustados
ggplot(data = datos2, aes(x = fitted(modelo_con_col), y = rstandard(modelo_con_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "valores ajustados", y = "residuos estandarizados",
       title = "Residuos estandarizados vs Valores predichos")

  • Residuos vs. cada variable regresora
ggplot(data = datos2, aes(x = x1, y = rstandard(modelo_con_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = expression(x[1]), y = "residuos estandarizados",
       title = expression(paste("Residuos estandarizados vs. ", x[1])))

ggplot(data = datos2, aes(x = x2, y = rstandard(modelo_con_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = expression(x[2]), y = "residuos estandarizados",
       title = expression(paste("Residuos estandarizados vs. ", x[2])))

ggplot(data = datos2, aes(x = x3, y = rstandard(modelo_con_col))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = expression(x[3]), y = "residuos estandarizados",
       title = expression(paste("Residuos estandarizados vs. ", x[3])))

15.9 Normalidad de los residuos

residuos.df=data.frame(resid=rstandard(modelo_con_col))
ggplot(data = residuos.df, aes(sample = rstandard(modelo_con_col))) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "QQ Plot de los residuos",
       x = "cuantiles teóricos", y = "cuantiles de los residuos") +
  theme_minimal()

16 Detectando multicolinealidad

library(car)

vif(modelo_sin_col)
##       x1       x2       x3 
## 1.019125 1.003057 1.017576
vif(modelo_con_col)
##         x1         x2         x3 
## 354.260148   1.023393 354.548013

16.1 ¿Cómo tratar el problema de multicolinealidad?

  • Lo más fácil es elegir un subconjunto de las variables regresoras poco correlacionadas.

  • El asunto es… si detectamos dos variables muy correlacionadas ¿cuál omito? En general, conviene omitir aquella variable que tenga:

    • muchos datos faltantes
    • mayor error de medición.
  • Otra posibilidad es utilizar métodos para eliminar variables a través de procedimientos de selección automáticos (se presentarán más adelante).

  • También se puede consturir una tercer variable que sea combinación de estas variables correlacionadas. Por ejemplo: el índice de precios.

16.2 Métodos para tratar el problema de la multicolinealidad

  • Selección de variables: con este método trabajaremos.

  • Métodos de regularización o penalización. Ridge o Lasso: mejoran la estabilidad de los coeficientes, especialmente cuando hay multicolinealidad o muchas variables. Añaden un término de penalización al criterio de ajuste del modelo. Es decir, en lugar de minimizar solo el error cuadrático (como en la regresión lineal clásica), también se penalizan los coeficientes grandes.

  • Métodos de reducción de la dimensión: son técnicas que reducen la cantidad de variables regresoras proyectando sobre un espacio de dimensión menor. Este método no se tratará en el curso.

17 Conclusión

Reformulamos el modelo eligiendo \(x_1\) o \(x_3\).


Referencias:

  • Montgomery, E., Vining, D. y Peck. (2006). Introducción Al Análisis de Regresión Lineal. 3ed ed. México: Cecsa.

  • James, G., Witten, D., Hastie, T., y Tibshirani, R. (2013). An introduction to statistical learning (Vol. 112). Springer.